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ABSTRACT 



The sample of known exoplanets is strongly biased to masses larger than 
the ones of the giant gaseous planets of the solar system. Recently, the 
discovery of two extrasolar planets of considerably lower masses around the 
nearby stars GJ 436 and p Cancri was reported. They are like our outer- 
most icy giants, Uranus and Neptune, but in contrast, these new planets 
are orbiting at only some hundredth of the Earth-Sun distance from their 
host stars, raising several new questions about their origin and constitution. 
Here we report numerical simulations of planetary accretion that show, for 
the first time through N-body integrations that the formation of compact 
systems of Neptune-like planets close to the hosts stars could be a common 
by-product of planetary formation. We found a regime of planetary accre- 
tion, in which orbital migration accumulates protoplanets in a narrow region 
around the inner edge of the nebula, where they collide each other giving rise 
to Neptune-like planets. Our results suggest that, if a protoplanetary solar 
environment is common in the galaxy, the discovery of a vast population of 
this sort of 'hot cores' should be expected in the near future. 
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1 Introduction 

The most plausible hypothesis regarding the formation of giant planets is 
furnished by the core accretion model (Mizuno, 1980), where a solid core 
grows by accretion of planetesimals (Safronov, 1969). When the core mass 
increases above a certain critical value (10-15 Earth masses), a violent ac- 
cretion of nebular gas onto the solid core starts, leading to the formation of 
a gaseous giant planet (Pollack et al, 1996). The final mass of the gaseous 
envelope depends on several factors, including the amount of gas present in 
the primordial disk, and the abihty of the solid core to reach the critical mass 
before the dissipation of the surrounding nebula. 

Two extrasolar planets of low mass were recently discovered orbiting 
around the nearby stars GJ 436 (Butler et al., 2004) and p Cancri (McArthur 
et ai, 2004). They are like Neptune in mass, but these new planets are or- 
biting at only 0.02 AU and 0.04 AU respectively from their host stars. 

Standard disk models show that at the distances where the new planets 
have been found (Butler et al., 2004; McArthur et al., 2004), the temperature 
is 600 - 2000 K (depending on the stellar type). Although at these distances, 
the existence of solid material cannot be ruled out, the temperature is too 
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hot as to expect a substantial amount of it as to form large solid cores there. 
More hkely, giant planets hke these ones form at much larger distances from 
the star, and subsequently they migrate inwards. Theory and numerical 
simulations have shown that growing protoplanets can experience orbital mi- 
gration, traveling very far from their birth place (Goldreich and Tremaine, 
1980; Lin and Papaloizou, 1986; Ward 1986; Ward 1997a). In fact, plan- 
etary migration is a natural consequence of the tidal interaction with the 
nebula. The most plausible mode of migration for protoplanets is inward in 
a quasi circular orbit (Ward, 1986; Ward, 1988; Ward, 1997a; Tanaka et al, 
2002; Tanaka and Ward, 2004). The tidal interaction with the star or the 
truncation of the disk by the stellar magnetosphere, prevent the planets to 
fall onto the star. Planetary migration have been invoked to explain close 
stellar companions (Lin et al, 1996) of large mass ('hot Jupiters'). Ward 
(1997b), have also discussed in situ formation of ice giants planets by migra- 
tion of embryos to the inner protoplanetary disk limit. Besides, Ida and Lin 
(2004) have delimited conditions for which ice giants could form in a variety 
of circumstellar disks. Laughlin et al. ( 2004a), show that whereas ice giants 
form readily around stars of all masses, Jupiter-like planets are difficult to 
form around the smallest stars. The new ice planets discovered support the 
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findings of Lauglilin et al. (2004a); indeed, GJ436 is a M dwarf type star, 
but in the p Cancri system (which is a 'solar' G8V star) are coexisting the 
new 'Neptune' with three other Jupiter-hke planets. However, N-body simu- 
lations of planetary formation showing how Neptune-like planets as the ones 
recently discovered could form, have not been yet carried out. In this paper, 
we report the first numerical simulations of planetary accretion, showing that 
the formation of Neptune-like planets very near the central star could be a 
common by product of planet formation, in a standard nebular environment. 

2 Simulations 

We have performed a series of N-body numerical simulations of the accretion 
of sohd cores, that include, for the first time, the dynamical effects of the 
gaseous environment on them, as predicted by the more recent results of the 
theory of planet-disk interaction (Ward, 1997a; Tanaka et al., 2002; Tanaka 
and Ward, 2004). In the set of simulations reported here, we have explored a 
model whose initial mass is of the order of the mass expected in the minimum- 
mass nebular model (NM) (Hayashi et al., 1985). 

For the accretion of solid cores, we follow the most plausible model which 
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is based on the concept of 'oligarchic' growth (Kokubo and Ida, 1998; Kokubo 
and Ida, 2000). In this model, only few large objects growth at almost compa- 
rable rates, separated by amounts determined by their masses and distances 
from the star. The gravitational scattering with the largest protoplanets 
dominates the dynamical evolution of the background planetesimals, which 
for this reason cease to growth. However, the planetesimal disk, although 
representing only a fraction of the total solid mass, may contribute to the 
further growth of migrating protoplanets. We have performed four numerical 
simulations where 100 small protoplanets of 0.5 Earth masses were initially 
placed on nearly circular, very low-inclination orbits around a star of one 
solar mass (M*). They were spread from 5 AU to 15 AU, their mutual sep- 
aration were generated at random, but following a r~^^^ profile. 

It is very difficult to handle the amount of solid material required for a 
NM (~ TOM®) -note that we must integrate bodies spread over an extended 
disk-; thus, we are compelled to make some approximations. Oligarchyc 
growth models predict (at an initial time) a variety of embryo masses at 
different heliocentric distances, but we assume a uni-modal embryonic mass 
distribution. The number and mass of the embryos cited above, are plausible 
initial conditions in an protoplanetary disk model of a few minimum-masses 

8 



(see, Thommes et ai, 2003). Although the time- forming of such embryos can 
be very large, even in a several NM model, if the random velocities of the 
embryos are sufficiently damped, the formation of big bodies can be expected 
(Thommes et ai, 2003). We make the hypothesis that the tidal interaction 
with a NM gaseous disk (not considered in ohgarchyc growth studies) can 
produce substantial velocity damping and radial mobility of big embrios as 
to strength the growth process (see e.g., Tanaka and Ida, 1999). Besides, 
during the epoch of the formation of these embryos the nebular environment 
could have been denser, helping to the growth process. Nevertheless, we 
think that detailed time and mass dependence of the initial conditions may 
not be important to a bottom-up set of simulations. 

In addition, a swarm of residual planetesimals was included, occupying 
the same space than the protoplanets. To make the simulations numerically 
tractable, we used 200 planetesimals of 0.1 Earth masses, distributed with 
a surface density proportional to r'^-^. Although the planetesimals interact 
with the protoplanets, we ignored the self-gravity between them (see e.g., 
Cionco and Brunini, 2002). A gaseous disk consistent with the minimum 
mass solar nebula model (total mass, density and temperature distribution, 
and scale height) was considered (Hayashi et ai, 1995). On the planetary 
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embryos, the disk acts dynamically on the orbital semimajor axes and eccen- 
tricities, as prescribed by the theory of density waves. The gaseous disk was 
assumed truncated at 0.1 AU and gradually dissipates at a constant rate, 
in such a way that it completely disappears after 10^ yr, as observations 
of circumstellar disks around young stars indicates (Beckwith and Sargent, 
1996). The gaseous disk also affects the orbits of the planetesimals through 
aerodynamic gas drag (Adachi et al., 1976), but as if their radii were of 1000 
km. 

To carry out these simulations, we have adapted our N-Body hybrid code 
(Brunini and Mehta, 2002; Cionco and Brunini, 2002) to include all these 
effects. Our model of planet-disk interaction includes migration of type I. 

We didn't follow the 'classical' approach of summing the corresponding 
component torque over each active resonance (see e.g., Cionco and Brunini, 
2002); here, we follow the most efficient strategy to be applied either in the 
leap-frog part of the hybrid integrator as in the Bulirsch-Stoer numerical 
scheme used in it, in order to prescribe the orbital evolution predicted by 
density waves theory. We correct the embryo velocities through a Stokes non- 
conservative force, yielding the same semi-major axis (a) and eccentricity (e) 
evolution predicted by density waves up to first order in e, applied in a 
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consistent fashion, taking into account the disk properties at the different 
zones reached by migrating bodies. 

This force, acting on a protoplanet of mass Mp and mean motion fl, 
produces an acceleration {As) and must be apphed following the gas-planet 
relative velocity direction (vrei) as in the case of aerodynamic drag (Adachi 
et ai, 1976): 

As^KMppc'^Vrei, (1) 

where K is a. parameter; p is the volumetric gaseous density; is the sound 
speed in the medium (assumed as gas temperature dependent); Vrei is the 
modulus of the relative gas-planet velocity. In what follows, we assume that 
all the distance-dependent quantities are evaluated at the semi- major axis of 
the perturber. 

Adachi et al. (1976), have found expressions for the mean variation of the 
orbital elements over one orbital period due to Eq. (1). Rewriting p oc 
being E the gaseous surface density, and h the isothermal vertical scale height 
defined by: h oc Cs/fl, we can write up to the first order in e: 
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U^\^-KM,T.ncr'{h/a)\ (2) 



-/^\ oc -i^MpEQcr^ (3) 
e \at I 

The perturbing acceleration -Eq. (1)-, must reproduce the same evolu- 
tion of semi-major axis and eccentricity estimated by density waves. This 
estimates are: (see e.g., Ward, 1986; Ward, 1988; Ward, 1997a; Tanaka et 
al., 2002; Tanaka and Ward, 2004): 



Ida Mp /Sa^\ ^ / a\ 

— r -ITT- -in- ^ > (4) 



edt ) ^ ' ' ^ ' 

Expressions (2)-(3) and (4)-(5) are equivalent adopting a — —3, and K oc 
G^, being G — Q?' /M^, the gravitational constant. 

Finally, we adopted an appropriated constant in the definition of the pa- 
rameter K in order to match 'exactly' the semi-major axis and eccentricity evolution 
result predicted by Eqs.(4)-(5) in our simulated disk. 
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In our simulations, a 1 planet migrating at 5 AU, reach the inner 
disk limit in ~ 1 x 10^ yr, and the time required for the same body, to reach 
a circular orbit (i.e., in attain e ~ 5 x 10~^), is of ~ 5 x 10^ yr. 

The inner edge of the disk is an uncertain parameter, nevertheless, the 
migration can be halted by gas depletion near 0.1 AU (see Ida and Lin, 2004, 
and references therein); thus, we assume that a protoplanet reaching the 
'real' disk limit, can go on migrating until its outer Lindblad resonances get 
out of the disk by a factor of ~ 2 respect to the physical disk limit. Then, 
according to the frictional approach, when a protoplanet reaches 0.05 AU, it 
is assumed that the inward migration stops. 

In these simulations, if a planetary core reaches 10 Earth masses, it start 
to accrete a gaseous envelope following the core-instability model (Pollack et 
al., 1996). The gas accretion is prescribed adding an additional mass to the 
embryo at each timestep, following the standard (Jl) model of Pollack et ai, 
1996. 

In the simulations when the numerical Bulirsch-Stoer part of the in- 
tegrator is used (i.e., to integrate accurately embryo-embryo or embryo- 
planetesimal close encounters, close stellar approaches, or the swarm of bod- 
ies reaching the inner disk limit), the time step is automatically adapted. 
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We use a tolerance of 10~^^. Other effects, as tidal interaction with the cen- 
tral star, quadrupole distortion between the planets, or relativistic effects, 
where not considered in these set of simulations (see e.g.. Trilling et al., 

1998, Kiseleva-Egletton and Bois, 2001). 

3 Results 

The combined results of the four simulations are shown in Fig. 1. It is evi- 
dent that several cores (only one of the planets accrete negligible amounts of 
gas) reach the inner edge of the simulated disk, where their inward migration 
stops, and then, they survive with the semi-major axis reported. The migra- 
tion of the cores is always of type I because their masses are not large enough 
as to open a gap in the disk. By 10 Myr, when the disk is assumed completely 
dispersed, several 'hot cores' resembling the ones found around GJ 436 and 
p Cancri are present. During the first stage of the simulations, they growth 
by the accretion of planetesimals, which is possible for two reasons: first, the 
fast type I orbital migration, which is shown in Fig. 2, allows the embryos to 
be continuously immersed in zones not previously depleted of planetesimals 
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^ and second, the aerodynamic drag maintains the planetesimals in relatively 
low excited orbits, a condition needed for efficient accretion. But the most 
important growth regime is due to mutual accretion between massive em- 
bryos during the last stages of the process. This happens because a fraction 
of them migrated quickly to a narrow zone around the inner edge of the 
disk favoring mutual collisions and the subsequent formation of Neptune-like 
planets. 

Table 1 shows the final 'hot cores' obteined in each simulation (only plan- 
ets with Mp > lOM® are included) at t = 10 Myr. The planets of the simula- 
tions Nr2 and Prl are remarkably similar to the ones found around GJ 436: 
21 Me with a = 0.0278, (Butler et a/., 2004); and p Cancri: 17.7 M® with 
a=0.04, (McArthur et al., 2004). In Fig. 3, we have reproduced the growing 
process of two planets in the above mentioned simulations. An interesting 
question arising from the simulations, is the amount of solid mass engulfed 
by the star. Although it is in general negligible, simulation Pr2, the only one 

with two final cores, disperses toward the star about one core mass of solid 
^Note that if the disk is turbulent, the cores can attain a random migration according 

to LaughUn et al. (2004b) not considered in our model; nevertheless, these fluctuations are 

largely contained in the random walk introduced by the perturbations between embryos. 
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material. We believe that our simulations are the first showing this growth 
regime, that operates on type I migrating protoplanets, in a standard neb- 
ular environment. If the new discovered planets formed in a scenario as the 
one described here, they should be of ice-rock composition with only a thin 
atmosphere. 

4 Conclusions 

Giant gaseous planets as massive as Jupiter, or even Saturn, were not found 
in this set of runs, but in all the four simulations, compact systems of two 
or more Neptune-like planets were produced. Although this result points to 
the existence of compact planetary systems formed by small giants, it is in 
contrast, however, to the system orbiting p Cancri, which already has other 
three Jupiter-like planets. This fact could be due to the imposed restriction 
on the critical core mass, for fast accretion of gas, which nevertheless is not 
the only possible scenario for the core instabihty model (Stevenson, 1982). 
In fact, according to recent results of hydrodynamic accretion of gas around 
solid cores (D'Angelo et al., 2002), planetary embryos may accrete gas when 
the solid cores reach only 1-2 Earth masses. 

16 



Our results show that the formation of planets, and systems of Neptune- 
like planets around solar type stars, as the ones recently discovered, seems to 
be a natural by-product of planetary formation, and encourage further stud- 
ies statistically more robust, exploring the whole set of parameters i nvolved 
in the problem. Planetary migration and the value of the critical core mass 
for substantial accretion of gas, are crucial factors in determining the distri- 
bution of orbits and masses of the formed planetary systems. Taking into 
account the standard scenario explored here, we predict that a substantial 
population of 'hot cores' with masses between ~ 14 and ~ 24 M® placed 
at semi-major axis< 0.05 AU await to be discovered in the near future. 
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referees that help us to improve the manuscript. 
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Figure captions: 

• Figure 1 Final planetary masses of the embryos as a function of the 
mean distance to the star. The different symbols represent the planets 
obtained in each simulation after 10 Myr. In the four simulations, 
planets of 10-24 Earth masses were produced very close to central star. 

• Figure 2 Evolution of the semimajor axis of one of the planetary cores 
in our simulations. The orbital migration obeys to the type 1 migra- 
tion, as the theory of disk planet interaction predicts. The migration 
almost stops when the core reaches the inner edge of the simulated 
disk. At this location the core suffers strong dynamical interactions 
with other embryos that also have reached the inner edge of the disk. 
These interactions produce the observed jumps after 5.5 My at short 
semimajor axes. 

• Figure 3 Mass accretion of two Neptune-like planets vs. the mean 
distance to the host star. Although the embryos accrete some back- 
ground planetesimals (small jumps in mass at large semimajor axes), 
the accretion of substantial amounts of mass is due to the mutual coag- 
ulation between embryos, represented by the large steps in mass once 
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they reach the inner edge of the disk. The presence of more than one 
embryo in a narrow zone around the inner edge of the disk produce the 
jumpy evolution in semimajor axis shown in Fig. 2. 
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Table 1: Results of the four simulations at 10 Myr: name of the simula- 
tion; total mass, gaseous mass, semi-major axis and eccentricity of the cores 
obtained; M^ng is the mass of solid material engulfed by the star. 
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Figure 1 
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